1 Loading the analytical data

overall <- read_rds(here("data", "output_data", "overall.rds"))

black_miss <- read_rds(here("data", "output_data", "black.rds")) 

black_complete <- read_rds(here("data", "output_data", "black.rds")) %>% 
  #removed states with partially missing data (i.e., between 4 and 19 yearly outcome data
  #missing)  #for four states (Alaska, New Mexico, Rhode Island, Utah) and completely
  #missing   #for nine states (Hawaii, Idaho, Maine, Montana, New Hampshire, North Dakota,
  #South Dakota, Vermont, Wyoming. Total removed (n=4 + 9 = 13)
  filter(!(state %in% c('Alaska', 'New Mexico','Rhode Island', 'Utah',
                        
                        'Hawaii', 'South Dakota',
                        'Idaho', 'Maine', 'Montana', 'New Hampshire',
                        'North Dakota', 'Vermont', 'Wyoming')))

hispanic_miss <- read_rds(here("data", "output_data", "hispanic.rds"))

hispanic_complete <- read_rds(here("data", "output_data", "hispanic.rds")) %>% 
  #removed states with partially missing data (i.e., between 4 and 19 yearly outcome data missing) 
  # for twelve (12) states (Alabama, Arkansas, Idaho, Iowa, Kentucky, Louisiana, 
  # Minnesota, Nebraska, Rhode Island, South Carolina, Tennessee, Wisconsin) and completely missing for 
  # eleven (11) states (Alaska, Delaware, Maine, Mississippi, Montana, New Hampshire, 
  # North Dakota, South Dakota, Vermont, West Virginia, Wyoming). Total removed (n=11+12 = 23)
  filter(!state %in% c("Alabama", "Arkansas", "Idaho", "Iowa", "Kentucky", 
                       "Louisiana", "Minnesota", "Nebraska", "Rhode Island",
                       "South Carolina", "Tennessee", "Wisconsin",
                       
                       "Alaska","Delaware","Maine","Mississippi","Montana",
                       "New Hampshire","North Dakota",  "South Dakota",  
                       "Vermont","West Virginia", "Wyoming")) %>% 
  #We imputed below outcome data for the six states with minimal missing (Kansas, 
  # Maryland, Missouri, North Carolina, Oregon, Utah). We imputed the 
  # outcome data with the closest outcome data (previous/next year) given 
  # that excluding these six states will have prevented the generalized 
  # synthetic control method from working
  arrange(state, year) %>% 
  dplyr::group_by(state) %>%
  fill(c("cvd_death_rate", "population"), .direction = "downup") %>%
  dplyr::ungroup() 

white <- read_rds(here("data", "output_data", "white.rds")) 

men <- read_rds(here("data", "output_data", "men.rds")) 

women <- read_rds(here("data", "output_data", "women.rds"))

2 Analyzing the data using the generalzied SCM

2.1 Overall

est_overall <- gsynth(cvd_death_rate ~ treatedpost + primarycare_rate + 
                        cardio_rate + population + 
                        low_educ + employed_for_wages + party +  
                        low_income + married + male + race_nonwhite,
                      data = overall,
                      EM = F, 
                      index = c("state_id","year"),
                      inference = "parametric", se = TRUE, 
                      r = c(0, 5), 
                      seed = 123,
                      nboots = 200, CV = TRUE, force = "two-way", parallel = FALSE)


#ATT plot
p1 <- plot(est_overall, type = "counterfactual", raw = "none", 
           theme.bw = TRUE, main = "", #ylim = c(120, 200), 
           legendOff = TRUE, xlab = "", ylab = "") +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        #axis.text.x = element_text(angle = 90, vjust = 0.5),
        plot.margin = margin(10, 10, 50, 10))
overall_by_period <- est_overall[["est.att"]] %>% 
  as.data.frame() %>% 
  rownames_to_column(var= "time") %>% 
  mutate(time=as.numeric(time)) %>% 
  dplyr::filter(time>=0)

overall_average <- data.frame(est_overall$est.avg)
overall_average

2.2 Black

est_black <- gsynth(cvd_death_rate ~ treatedpost + primarycare_rate + 
                        cardio_rate + population + 
                        low_educ + employed_for_wages + party +  
                        low_income + married + male,
                    data = black_complete, 
                    EM = F, 
                    index = c("state_id","year"),
                    inference = "parametric", se = TRUE, 
                    r = c(0, 5), 
                    seed = 123,
                    nboots = 200, CV = TRUE, force = "two-way", parallel = FALSE)


p2 <- plot(est_black, type = "counterfactual", raw = "none", theme.bw = TRUE, main = "",#ylim = c(220, 475),  
           legendOff = TRUE, xlab = "", ylab = "") +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        #axis.text.x = element_text(angle = 90, vjust = 0.5),
        plot.margin = margin(10, 10, 50, 10))
black_by_period <- est_black[["est.att"]] %>% 
  as.data.frame() %>% 
  rownames_to_column(var= "time") %>% 
  mutate(time=as.numeric(time)) %>% 
  dplyr::filter(time>=0)

black_average <- data.frame(est_black$est.avg)
black_average

2.3 Hispanic

est_hispanic <- gsynth(cvd_death_rate ~ treatedpost + primarycare_rate + 
                        cardio_rate + population + 
                        low_educ + employed_for_wages + party +  
                        low_income + married + male,
                       data = hispanic_complete,
                       EM = F,
                       index = c("state_id","year"),
                       inference = "parametric", se = TRUE,
                       r = c(0, 5),
                       seed = 123,
                       nboots = 200, CV = TRUE, force = "two-way", parallel = FALSE)

p3 <- plot(est_hispanic, type = "counterfactual", raw = "none", theme.bw = TRUE, main = "",#ylim = c(220, 475),
           legendOff = TRUE, xlab = "", ylab = "") +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        #axis.text.x = element_text(angle = 90, vjust = 0.5),
        plot.margin = margin(10, 10, 50, 10))
hispanic_by_period <- est_hispanic[["est.att"]] %>%
  as.data.frame() %>%
  rownames_to_column(var= "time") %>%
  mutate(time=as.numeric(time)) %>%
  dplyr::filter(time>=0)

hispanic_average <- data.frame(est_hispanic$est.avg)
hispanic_average

2.4 White

est_white <- gsynth(cvd_death_rate ~ treatedpost + primarycare_rate + 
                        cardio_rate + population + 
                        low_educ + employed_for_wages + party +  
                        low_income + married + male,
                    data = white,  
                    EM = F, 
                    index = c("state_id","year"),
                    inference = "parametric", se = TRUE, 
                    r = c(0, 5), 
                    seed = 123,
                    nboots = 200, CV = TRUE, force = "two-way", parallel = FALSE)

p4 <- plot(est_white, type = "counterfactual", raw = "none", theme.bw = TRUE, main = "",#ylim = c(220, 475),  
           legendOff = TRUE, xlab = "", ylab = "") +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        #axis.text.x = element_text(angle = 90, vjust = 0.5),
        plot.margin = margin(10, 10, 50, 10))
white_by_period <- est_white[["est.att"]] %>% 
  as.data.frame() %>% 
  rownames_to_column(var= "time") %>% 
  mutate(time=as.numeric(time)) %>% 
  dplyr::filter(time>=0)

white_average <- data.frame(est_white$est.avg)
white_average

2.5 Men

est_men <- gsynth(cvd_death_rate ~ treatedpost + primarycare_rate + 
                        cardio_rate + population + 
                        low_educ + employed_for_wages + party +  
                        low_income + married + race_nonwhite,
                  data = men,  
                  EM = F, 
                  index = c("state_id","year"),
                  inference = "parametric", se = TRUE, 
                  r = c(0, 5), 
                  seed = 123,
                  nboots = 200, CV = TRUE, force = "two-way", parallel = FALSE)

p5 <- plot(est_men, type = "counterfactual", raw = "none", theme.bw = TRUE, main = "",#ylim = c(220, 475),  
           legendOff = TRUE, xlab = "", ylab = "") +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        #axis.text.x = element_text(angle = 90, vjust = 0.5),
        plot.margin = margin(10, 10, 50, 10))
men_by_period <- est_men[["est.att"]] %>% 
  as.data.frame() %>% 
  rownames_to_column(var= "time") %>% 
  mutate(time=as.numeric(time)) %>% 
  dplyr::filter(time>=0)

men_average <- data.frame(est_men$est.avg)
men_average

2.6 Women

est_women <- gsynth(cvd_death_rate ~ treatedpost + primarycare_rate + 
                        cardio_rate + population + 
                        low_educ + employed_for_wages + party +  
                        low_income + married + race_nonwhite,
                    data = women,  
                    EM = F, 
                    index = c("state_id","year"),
                    inference = "parametric", se = TRUE, 
                    r = c(0, 5), 
                    seed = 123,
                    nboots = 200, CV = TRUE, force = "two-way", parallel = FALSE)

p6 <- plot(est_women, type = "counterfactual", raw = "none", theme.bw = TRUE, main = "",#ylim = c(220, 475),  
           legendOff = TRUE, xlab = "", ylab = "") +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        #axis.text.x = element_text(angle = 90, vjust = 0.5),
        plot.margin = margin(10, 10, 50, 10))
women_by_period <- est_women[["est.att"]] %>% 
  as.data.frame() %>% 
  rownames_to_column(var= "time") %>% 
  mutate(time=as.numeric(time)) %>% 
  dplyr::filter(time>=0)

women_average <- data.frame(est_women$est.avg)
women_average

3 Manuscript Figures and Tables

3.1 Table 1: Baseline characteristics

Baseline characheristics

## Table 1: Baseline Characteristics----
#overall <- read_rds(here("data", "output_data", "overall.rds"))
tab <- overall %>% 
  filter(year==2014) %>% 
  set_variable_labels(
    cvd_death_rate = "CVD deaths per 100,000 persons aged 45-64 years",
    primarycare_rate = "Primary care clinicians per 100,000 residents",
    cardio_rate = "Cardiologists per 100,000 residents",
    low_income = "Percentage of residents aged 45-64 years with annual income less than $15,000",
    male = "Percentage of residents aged 45-64 years who are males",
    race_nonwhite = "Percentage of residents aged 45-64 years who are Non-White",
    married = "Percentage of residents aged 45-64 years who are married",
    low_educ = "Percentage of residents aged 45-64 years without high school degree",
    employed_for_wages = "Percentage of residents aged 45-64 years who are employed for wages",
    party = "Percentage of political party") %>% 
  select(cvd_death_rate, primarycare_rate, cardio_rate,low_income,
         male, race_nonwhite, married, low_educ, employed_for_wages,
         party, treated) %>% 
  mutate_at(vars(low_educ, employed_for_wages, 
                 low_income, married, male, race_nonwhite), multiply_by, e2=100) %>% 
  set_value_labels(treated = c("Medicaid Non-Expansion States"= 0,
                               "Medicaid Expansion States"=1),
                   party = c("Republican"=0, "Democrat"=1, "Split"=2)) %>% 
  modify_if(is.labelled, to_factor) %>% 
  tbl_summary(by = treated,
              missing = 'no',
              digits = list(all_continuous()~1,
                            all_categorical()~1), 
              statistic = list(all_continuous() ~ '{mean} ({sd})',
                               all_categorical() ~ '{n} ({p}%)')) %>% 
    modify_header(label = "**Characteristics**") %>%
    add_overall() %>%
    add_stat_label(location = "row") %>%
    modify_spanning_header(c("stat_1","stat_2") ~ "**State Medicaid Expansion Status**") %>% bold_labels()

tab
Characteristics Overall, N = 50 State Medicaid Expansion Status
Medicaid Non-Expansion States, N = 16 Medicaid Expansion States, N = 34
CVD deaths per 100,000 persons aged 45-64 years, Mean (SD) 154.0 (47.0) 177.0 (52.5) 143.1 (40.5)
Primary care clinicians per 100,000 residents, Mean (SD) 77.0 (12.8) 67.6 (7.2) 81.5 (12.5)
Cardiologists per 100,000 residents, Mean (SD) 6.1 (2.2) 5.2 (1.5) 6.5 (2.3)
Percentage of residents aged 45-64 years with annual income less than $15,000, Mean (SD) 9.6 (3.1) 10.7 (3.3) 9.1 (2.8)
Percentage of residents aged 45-64 years who are males, Mean (SD) 42.4 (2.9) 41.8 (3.1) 42.7 (2.8)
Percentage of residents aged 45-64 years who are Non-White, Mean (SD) 21.4 (12.9) 23.2 (12.5) 20.5 (13.3)
Percentage of residents aged 45-64 years who are married, Mean (SD) 61.4 (4.4) 61.9 (5.6) 61.2 (3.8)
Percentage of residents aged 45-64 years without high school degree, Mean (SD) 6.8 (3.1) 8.0 (3.5) 6.2 (2.8)
Percentage of residents aged 45-64 years who are employed for wages, Mean (SD) 63.9 (6.9) 61.8 (7.8) 64.9 (6.2)
Percentage of political party, n (%)
    Republican 24.0 (48.0%) 15.0 (93.8%) 9.0 (26.5%)
    Democrat 13.0 (26.0%) 0.0 (0.0%) 13.0 (38.2%)
    Split 13.0 (26.0%) 1.0 (6.2%) 12.0 (35.3%)

3.2 Table 2: Overall and subgroup effect

Overall and subgroup effect of the Medicaid expansion on CVD mortality

overall_average$group='overall'

black_average$group='black'

hispanic_average$group='hispanic'

white_average$group='white'

men_average$group='men'

women_average$group='women'


average_data <- rbind(overall_average, black_average, hispanic_average, 
                 white_average, men_average, women_average)

average_data$group <- factor(average_data$group, 
                        levels=c('overall', 'black', 'hispanic', 'white', 
                                 'men','women'))

average_table <- average_data %>% 
  transmute(Group = group,
            `Adjusted Mean Difference (95%CI)`= paste(digits(Estimate,2),"", " " ,"(",
                                                   digits(CI.lower,2),", ",
                                                   digits(CI.upper,2), ")", sep=""))
            

knitr::kable(average_table, caption = "Overall and subgroup effect of the Medicaid expansion on CVD mortality") %>%
  kableExtra::kable_styling(bootstrap_options = c("striped", "bordered"),
                full_width = F)
Overall and subgroup effect of the Medicaid expansion on CVD mortality
Group Adjusted Mean Difference (95%CI)
ATT.avg overall -4.29 (-9.87, 1.29)
ATT.avg1 black -5.36 (-22.63, 11.91)
ATT.avg2 hispanic -4.28 (-30.08, 21.52)
ATT.avg3 white -3.18 (-8.30, 1.94)
ATT.avg4 men -5.96 (-15.42, 3.50)
ATT.avg5 women -3.34 (-8.05, 1.37)
# tripple difference in difference

ddd_black_hispanic_vs_white_overall <- bind_rows(black_average,
                                         hispanic_average,
                                         white_average) %>% 
  dplyr::select(group, ATT=Estimate, CI.lower, CI.upper, se=S.E.) %>% 
  pivot_wider(., names_from=group, 
              values_from=c(ATT, CI.lower, CI.upper, se)) %>% 
  summarise(mean_black=round((ATT_black-ATT_white), 2),
            se_black=sqrt(se_black^2 + se_white^2),
            lower_black=round((mean_black-qnorm(p=1-(0.05)/2, mean=0, sd=1)*se_black), 2),
            upper_black=round((mean_black+qnorm(p=1-(0.05)/2, mean=0, sd=1)*se_black), 2),
            mean_hispanic=round((ATT_hispanic-ATT_white), 2),
            se_hispanic=sqrt(se_hispanic^2 + se_white^2),
            lower_hispanic=round((mean_hispanic-qnorm(p=1-(0.05)/2, mean=0, sd=1)*se_hispanic), 2),
            upper_hispanic=round((mean_hispanic+qnorm(p=1-(0.05)/2, mean=0, sd=1)*se_hispanic), 2)) %>% 
   pivot_longer(cols=everything(),
                names_to = c(".value", "race"),
                names_sep="_") %>% 
  arrange(race)



ddd_women_vs_men_overall <- bind_rows(men_average,
                                        women_average) %>% 
  dplyr::select(group, ATT=Estimate, CI.lower, CI.upper, se=S.E.) %>% 
  pivot_wider(., names_from=group, 
              values_from=c(ATT, CI.lower, CI.upper, se)) %>% 
  summarise(mean=round((ATT_women-ATT_men), 2),
            se=sqrt(se_men^2 + se_women^2),
            lower=round((mean-qnorm(p=1-(0.05)/2, mean=0, sd=1)*se), 2),
            upper=round((mean+qnorm(p=1-(0.05)/2, mean=0, sd=1)*se), 2))


bind_rows(ddd_black_hispanic_vs_white_overall, ddd_women_vs_men_overall) %>% 
  rename(group = race) %>% 
  mutate(group = case_when(group=="black"~"Black vs White",
                           group=="hispanic"~"Hispanic vs White",
                           TRUE~"Women vs men")) %>% 
  mutate(`Adjusted Difference in Mean Difference (95%CI)`= paste(digits(mean,2),"", " " ,"(",
                                                                 digits(lower,2),", ",
                                                                 digits(upper,2), ")", sep="")) %>% 
  dplyr::select(Group= group,`Adjusted Difference in Mean Difference (95%CI)`) %>% 
  knitr::kable(caption = "Difference in mean differences") %>%
  kableExtra::kable_styling(bootstrap_options = c("striped", "bordered"),
                full_width = F)
Difference in mean differences
Group Adjusted Difference in Mean Difference (95%CI)
Black vs White -2.18 (-20.19, 15.83)
Hispanic vs White -1.10 (-27.40, 25.20)
Women vs men 2.62 (-7.95, 13.19)

3.3 Figure 1: Analytical data structure

Analytical data structure of the Medicaid expansion states and control states for the overall study population.

q1 <- panelview(cvd_death_rate ~treatedpost , data=overall, index=c("state","year"), 
                pre.post= TRUE,
                cex.axis.x=20,
                cex.axis.y=20,
                cex.lab=20)

3.4 Figure 2: Annual effect

Overall annual effect of Medicaid expansion on CVD deaths per 100,000 persons overall and by race and sex.

f_overall <- overall_by_period %>% 
  ggplot(aes(x = time, y = ATT)) +
  geom_point(color="#030303") + 
  geom_errorbar(aes( ymax = CI.lower, ymin = CI.upper) , width = .2, size = 0.7, 
                color="#030303")+
  geom_line(color="#030303") +
  ylim(-38, 38)+
  scale_x_continuous(breaks=seq(0,6,1))+
  geom_hline(yintercept = 0, lty=2,color="#030303") +
  ggtitle('Overall')+
  labs(x = "Number of years since Medicaid Expansion",
       y = "Change in CVD mortality, per 100,000 population") +
  theme_classic(base_size = 15) +
  theme(axis.line = element_line(colour = "grey50", size = 1),
        panel.background = element_blank(),
        axis.title = element_blank(),
        plot.title = element_text(hjust = 0.5))


f_black_hispanic_white <- bind_rows(black_by_period %>% mutate(group="black"),
                                    hispanic_by_period %>% mutate(group="hispanic"),
                                    white_by_period %>% mutate(group="white")) %>% 
  ggplot(aes(x = time, y = ATT, color = group)) +
  scale_color_manual(values=c("#8A2BE2", "#FF7F50", "#458B00"), 
                     labels=c('Black',   'Hispanic','White'))+
  geom_point(position = position_dodge(width = .5)) + 
  geom_errorbar( aes( ymax = CI.upper, ymin = CI.lower) , width = .2, 
                 position = position_dodge(width = .5), size = 0.7 )+
  geom_line(position = position_dodge(width = .5), size = 0.7) +
  ylim(-38, 38)+
  scale_x_continuous(breaks=seq(0,6,1))+
  #ggtitle('B')+
  geom_hline(yintercept = 0, lty=2, color="black") +
  theme(legend.position = "top",
        axis.title = element_blank(),
        axis.line = element_line(colour = "grey50", size = 1),
        panel.background = element_blank(),
        legend.background = element_rect(fill = "lemonchiffon", 
                                         colour = "grey50", 
                                         size = 1),
        legend.title = element_blank())

f_men_women <- bind_rows(men_by_period %>% mutate(group="men"),
                         women_by_period %>% mutate(group="women")) %>% 
  ggplot(aes(x = time, y = ATT, color = group)) +
  scale_color_manual(values=c("#1874CD", "#EE3B3B"), 
                     labels=c('Men','Women'))+
  geom_point(position = position_dodge(width = .5)) + 
  geom_errorbar( aes( ymax = CI.upper, ymin = CI.lower) , width = .2, 
                 position = position_dodge(width = .5), size = 0.7 )+
  geom_line(position = position_dodge(width = .5), size = 0.7) +
  ylim(-38, 38)+
  scale_x_continuous(breaks=seq(0,6,1))+
  #ggtitle('B')+
  geom_hline(yintercept = 0, lty=2, color="black") +
  theme(legend.position = "top",
        axis.title = element_blank(),
        axis.line = element_line(colour = "grey50", size = 1),
        panel.background = element_blank(),
        legend.background = element_rect(fill = "lemonchiffon", 
                                         colour = "grey50", 
                                         size = 1),
        legend.title = element_blank())

f_by_overall_race_sex <- ggarrange(f_overall, f_black_hispanic_white, f_men_women,
                                   ncol = 1, nrow=3) %>% 
  annotate_figure(left = textGrob("Change in CVD mortality, per 100,000 population", 
                                  rot = 90, vjust = 1, gp = gpar(cex = 1.3)), #cex = 0.8
                  bottom = textGrob("Number of years since Medicaid Expansion", gp = gpar(cex = 1.3)))

f_by_overall_race_sex

## Figure 3. Pre-treatment fit
Pre-treatment fit showing the observed CVD deaths per 100,000 persons along with the counterfactual (synthetic control). Vertical line represents the beginning of expansion. The solid line is expansion state, dashed line is the synthetic control.

p1_6_fit <- ggarrange(
  p1+ rremove("ylab") + rremove("xlab")+ggtitle("Overall"), 
  p2+ rremove("ylab") + rremove("xlab")+ggtitle("Black"),
  p3+ rremove("ylab") + rremove("xlab")+ggtitle("Hispanic"),
  p4+ rremove("ylab") + rremove("xlab")+ggtitle("White"),
  p5+ rremove("ylab") + rremove("xlab")+ggtitle("Men"),
  p6+ rremove("ylab") + rremove("xlab")+ggtitle("Women"),
  font.label = list(size = 2, color = "black", face = "bold", family = NULL),
  ncol=2, nrow=3,common.legend = TRUE)


p1_6_fit

3.5 eTable 1: STROBE Statement—Checklist (See Appendix)

3.6 eTable 2: Numeric values for Figure 1

Overall annual effect of the Medicaid expansion on CVD deaths per 100,000 persons for overall population and for the Black, Hispanic, White, Men and women subpopulations.

average_by_period_table <- 
  bind_rows(overall_by_period %>%  mutate(group="Overall"),
            black_by_period %>% mutate(group="Black"),
            hispanic_by_period %>% mutate(group="Hispanic"),
            white_by_period %>% mutate(group="White"),
            men_by_period %>% mutate(group="Men"),
            women_by_period %>% mutate(group="Women")) %>% 
  mutate(`Adjusted Mean Difference (95%CI)`= paste(digits(ATT,2),"", " " ,"(",
                                                   digits(CI.lower,2),", ",
                                                   digits(CI.upper,2), ")", sep="")) %>% 
  dplyr::select(group, time, `Adjusted Mean Difference (95%CI)`)
average_by_period_table %>% 
  knitr::kable(caption = "Overall annual effect of the Medicaid expansion on CVD deaths per 100,000 persons for the overall population and for the Black, Hispanic, White, Men and women subpopulations") %>%
  kableExtra::kable_styling(bootstrap_options = c("striped", "bordered"),
                full_width = F)
Overall annual effect of the Medicaid expansion on CVD deaths per 100,000 persons for the overall population and for the Black, Hispanic, White, Men and women subpopulations
group time Adjusted Mean Difference (95%CI)
Overall 0 0.40 (-3.38, 4.17)
Overall 1 -5.09 (-9.34, -0.85)
Overall 2 -7.61 (-13.43, -1.79)
Overall 3 -2.81 (-9.33, 3.71)
Overall 4 -4.18 (-10.76, 2.41)
Overall 5 -3.98 (-11.23, 3.28)
Overall 6 -1.55 (-9.88, 6.77)
Black 0 -1.48 (-16.15, 13.19)
Black 1 -6.19 (-23.10, 10.72)
Black 2 -14.67 (-32.74, 3.40)
Black 3 -4.66 (-23.60, 14.29)
Black 4 -5.44 (-26.24, 15.37)
Black 5 -1.61 (-25.21, 22.00)
Black 6 1.45 (-24.93, 27.82)
Hispanic 0 -2.56 (-18.05, 12.93)
Hispanic 1 -6.31 (-27.21, 14.60)
Hispanic 2 -7.93 (-31.09, 15.23)
Hispanic 3 -4.26 (-31.20, 22.68)
Hispanic 4 0.32 (-28.57, 29.21)
Hispanic 5 -6.78 (-37.26, 23.70)
Hispanic 6 -0.18 (-34.13, 33.78)
White 0 0.02 (-3.63, 3.68)
White 1 -5.42 (-10.99, 0.15)
White 2 -4.48 (-10.09, 1.13)
White 3 -2.16 (-8.10, 3.77)
White 4 -1.88 (-7.63, 3.87)
White 5 -3.06 (-9.81, 3.69)
White 6 -1.67 (-9.49, 6.15)
Men 0 1.50 (-5.36, 8.36)
Men 1 -4.72 (-12.59, 3.15)
Men 2 -10.35 (-20.28, -0.42)
Men 3 -2.40 (-12.92, 8.13)
Men 4 -8.16 (-19.81, 3.49)
Men 5 -5.74 (-18.13, 6.65)
Men 6 -4.16 (-17.59, 9.28)
Women 0 -1.58 (-6.37, 3.20)
Women 1 -5.50 (-9.92, -1.09)
Women 2 -4.36 (-9.87, 1.15)
Women 3 -2.91 (-8.15, 2.34)
Women 4 -2.06 (-7.69, 3.58)
Women 5 -4.57 (-10.57, 1.44)
Women 6 -0.08 (-5.49, 5.34)

3.7 eTable 3: Sensitivity analysis

Sensitivity analysis evaluating the effect of using heterogenous samples vs homogenous samples on both adjusted mean differences and difference in mean differences.

#Whites with black sample (n=37 states included)
white_withblacksample <- read_rds(here("data", "output_data", "white.rds")) %>% 
    filter(!(state %in% c('Alaska', 'New Mexico','Rhode Island', 'Utah',
                        
                        'Hawaii', 'South Dakota',
                        'Idaho', 'Maine', 'Montana', 'New Hampshire',
                        'North Dakota', 'Vermont', 'Wyoming')))

est_white_withblacksample <- gsynth(cvd_death_rate ~ treatedpost + primarycare_rate + 
                         cardio_rate + population + 
                         low_educ + employed_for_wages + party +  
                         low_income + married + male,
                       data = white_withblacksample,
                       EM = F,
                       index = c("state_id","year"),
                       inference = "parametric", se = TRUE,
                       #min.T0 = 5,
                       r = c(0, 5),
                       seed = 123,
                       nboots = 200, CV = TRUE, force = "two-way", parallel = FALSE)
## Cross-validating ... 
##  r = 0; sigma2 = 41.47336; IC = 4.63849; PC = 34.56113; MSPE = 52.04339
##  r = 1; sigma2 = 28.19493; IC = 4.91483; PC = 133.70899; MSPE = 37.70570*
##  r = 2; sigma2 = 21.85741; IC = 5.27680; PC = 189.27671; MSPE = 40.22883
##  r = 3; sigma2 = 17.07014; IC = 5.60049; PC = 214.83195; MSPE = 43.51960
##  r = 4; sigma2 = 15.45230; IC = 6.02614; PC = 255.26004; MSPE = 49.83340
##  r = 5; sigma2 = 13.26928; IC = 6.35339; PC = 271.50968; MSPE = 74.10160
## 
##  r* = 1
## 
## 
Simulating errors .....
Bootstrapping ...
## ..
est_white_withblacksample_average <- data.frame(est_white_withblacksample$est.avg)

est_white_withblacksample_average %>% 
  knitr::kable(caption = "Sensitivity analysis for the adjusted mean difference: Whites and Blacks have same sample (n=37)") %>%
  kableExtra::kable_styling(bootstrap_options = c("striped", "bordered"),
                full_width = F)
Sensitivity analysis for the adjusted mean difference: Whites and Blacks have same sample (n=37)
Estimate S.E. CI.lower CI.upper p.value
ATT.avg -2.15793 3.740264 -9.488713 5.172853 0.563976
#Whites with Hispanic sample (n=27 states included)

white_withhispanicsample <- read_rds(here("data", "output_data", "white.rds")) %>% 
  filter(!state %in% c("Alabama", "Arkansas", "Idaho", "Iowa", "Kentucky", 
                       "Louisiana", "Minnesota", "Nebraska", "Rhode Island",
                       "South Carolina", "Tennessee", "Wisconsin",
                       
                       "Alaska","Delaware","Maine","Mississippi","Montana",
                       "New Hampshire","North Dakota",  "South Dakota",  
                       "Vermont","West Virginia", "Wyoming")) 

est_white_withhispanicsample<- gsynth(cvd_death_rate ~ treatedpost + primarycare_rate + 
                                      cardio_rate + population + 
                                      low_educ + employed_for_wages + party +  
                                      low_income + married + male,
                                    data = white_withhispanicsample,
                                    EM = F,
                                    index = c("state_id","year"),
                                    inference = "parametric", se = TRUE,
                                    #min.T0 = 5,
                                    r = c(0, 5),
                                    seed = 123,
                                    nboots = 200, CV = TRUE, force = "two-way", parallel = FALSE)
## Cross-validating ... 
##  r = 0; sigma2 = 27.87331; IC = 4.56308; PC = 20.90499; MSPE = 55.58006*
##  r = 1; sigma2 = 16.13953; IC = 4.86382; PC = 141.06632; MSPE = 60.82607
##  r = 2; sigma2 = 12.59613; IC = 5.39248; PC = 210.92380; MSPE = 63.71601
##  r = 3; sigma2 = 7.25834; IC = 5.54719; PC = 179.74642; MSPE = 73.55888
##  r = 4; sigma2 = 4.99104; IC = 5.80804; PC = 163.69302; MSPE = 76.60868
##  r = 5; sigma2 = 2.24814; IC = 5.57526; PC = 91.82514; MSPE = 122.77355
## 
##  r* = 0
## 
## 
Simulating errors .....
Bootstrapping ...
## ..
est_white_withhispanicsample_average <- data.frame(est_white_withhispanicsample$est.avg)

est_white_withhispanicsample_average %>% 
  knitr::kable(caption = "Sensitivity analysis for the adjusted mean difference: Whites and Hispanics have same sample (n=27)") %>%
  kableExtra::kable_styling(bootstrap_options = c("striped", "bordered"),
                full_width = F)
Sensitivity analysis for the adjusted mean difference: Whites and Hispanics have same sample (n=27)
Estimate S.E. CI.lower CI.upper p.value
ATT.avg -11.44491 3.651427 -18.60158 -4.288245 0.0017223
ddd_black_vs_white_overall_homosample <- bind_rows(black_average %>%
                                                     mutate(group="black"),
  est_white_withblacksample_average %>%
  mutate(group="white")) %>% 
  dplyr::select(group, ATT=Estimate, CI.lower, CI.upper, se=S.E.) %>% 
  pivot_wider(., names_from=group, 
              values_from=c(ATT, CI.lower, CI.upper, se)) %>% 
  summarise(mean=round((ATT_black-ATT_white), 2),
            se=sqrt(se_black^2 + se_white^2),
            lower=round((mean-qnorm(p=1-(0.05)/2, mean=0, sd=1)*se), 2),
            upper=round((mean+qnorm(p=1-(0.05)/2, mean=0, sd=1)*se), 2),
            pvalue=round(pnorm(0,mean=abs(mean),sd=se) + 
                           (1 - pnorm(0,mean=-abs(mean),sd=se)),2))

ddd_black_vs_white_overall_homosample %>% 
  knitr::kable(caption = "Sensitivity analysis for the difference in mean differences: Whites and Blacks have same sample (n=37)") %>%
  kableExtra::kable_styling(bootstrap_options = c("striped", "bordered"),
                full_width = F)
Sensitivity analysis for the difference in mean differences: Whites and Blacks have same sample (n=37)
mean se lower upper pvalue
-3.2 9.570525 -21.96 15.56 0.74
ddd_hispanic_vs_white_overall_homosample <- bind_rows(hispanic_average %>%
                                                        mutate(group="hispanic"),
  est_white_withhispanicsample_average %>% mutate(group="white")) %>% 
  dplyr::select(group, ATT=Estimate, CI.lower, CI.upper, se=S.E.) %>% 
  pivot_wider(., names_from=group, 
              values_from=c(ATT, CI.lower, CI.upper, se)) %>% 
  summarise(mean=round((ATT_hispanic-ATT_white), 2),
            se=sqrt(se_hispanic^2 + se_white^2),
            lower=round((mean-qnorm(p=1-(0.05)/2, mean=0, sd=1)*se), 2),
            upper=round((mean+qnorm(p=1-(0.05)/2, mean=0, sd=1)*se), 2),
            pvalue=round(pnorm(0,mean=abs(mean),sd=se) + 
                           (1 - pnorm(0,mean=-abs(mean),sd=se)),2))

ddd_hispanic_vs_white_overall_homosample %>% 
  knitr::kable(caption = "Sensitivity analysis for the difference in mean differences: Whites and Hispanics have same sample (n=27)") %>%
  kableExtra::kable_styling(bootstrap_options = c("striped", "bordered"),
                full_width = F)
Sensitivity analysis for the difference in mean differences: Whites and Hispanics have same sample (n=27)
mean se lower upper pvalue
7.17 13.65951 -19.6 33.94 0.6

3.8 eFigure 1: Missing data Data structure for Blacks

Analytical data structure of Medicaid expansion states and control states for Blacks (before and after removing states with missing outcomes.

qblack_miss <- panelview(cvd_death_rate ~treatedpost , data=black_miss, index=c("state","year"), pre.post= TRUE)

qblack_complete <- panelview(cvd_death_rate ~treatedpost , data=black_complete, index=c("state","year"), pre.post= TRUE)

ggarrange(qblack_miss+ rremove("ylab") + rremove("xlab")+ggtitle("Black (with missing data)"), 
          qblack_complete+ rremove("ylab") + rremove("xlab")+ggtitle("Black (with complete data)"),
          ncol=1, nrow=2,
          common.legend = TRUE)

3.9 eFigure 2: Missing data Data structure for Hispanics

Analytical data structure of Medicaid expansion states and control states for Hispanics (before and after removing states with missing outcomes and nearest year imputation).

qhispanic_miss <- panelview(cvd_death_rate ~treatedpost , 
                            data=hispanic_miss, index=c("state","year"), 
                            pre.post= TRUE)

qhispanic_complete <- panelview(cvd_death_rate ~treatedpost , data=hispanic_complete, index=c("state","year"), pre.post= TRUE)

ggarrange(qhispanic_miss+ rremove("ylab") + rremove("xlab")+ggtitle("Hispanic (with missing data)"), 
          qhispanic_complete+ rremove("ylab") + rremove("xlab")+ggtitle("Hispanic (with complete and imputed data)"),
          ncol=1, nrow=2,
          common.legend = TRUE)

3.10 eFigure 3:Annual triple difference over time

Annual difference in mean difference between the effect of the Medicaid expansion on CVD deaths per 100,000 persons

## Black or Hispanic vs White-

ddd_black_hispanic_vs_white_by_period <- bind_rows(black_by_period %>% mutate(group="black"),
                hispanic_by_period %>% mutate(group="hispanic"),
                white_by_period %>% mutate(group="white")) %>% 
  dplyr::select(time, group, ATT, CI.lower, CI.upper, se=S.E.) %>% 
  pivot_wider(., names_from=group, 
              values_from=c(ATT, CI.lower, CI.upper, se)) %>% 
  summarise(mean_black=round((ATT_black-ATT_white), 2),
            se_black=sqrt(se_black^2 + se_white^2),
            lower_black=round((mean_black-qnorm(p=1-(0.05)/2, mean=0, sd=1)*se_black), 2),
            upper_black=round((mean_black+qnorm(p=1-(0.05)/2, mean=0, sd=1)*se_black), 2),
            mean_hispanic=round((ATT_hispanic-ATT_white), 2),
            se_hispanic=sqrt(se_hispanic^2 + se_white^2),
            lower_hispanic=round((mean_hispanic-qnorm(p=1-(0.05)/2, mean=0, sd=1)*se_hispanic), 2),
            upper_hispanic=round((mean_hispanic+qnorm(p=1-(0.05)/2, mean=0, sd=1)*se_hispanic), 2),
            time=time) %>% 
  relocate(time) %>% pivot_longer(cols=-time,
               names_to = c(".value", "race"),
               names_sep="_") %>% 
  arrange(race)

s1 <- ddd_black_hispanic_vs_white_by_period %>% 
  ggplot(aes(x = time, y = mean, color = race)) +
  scale_color_manual(values=c('#e41a1c', '#377eb8'), labels=c('Black vs White', 'Hispanic vs White'))+
  geom_point(position = position_dodge(width = .5)) + 
  geom_errorbar( aes( ymax = upper, ymin = lower) , width = .2, position = position_dodge(width = .5), size = 0.7 )+
  geom_line(position = position_dodge(width = .5), size = 0.7) +
  ylim(-38, 38)+
  scale_x_continuous(breaks=seq(0,6,1))+
  geom_hline(yintercept = 0, lty=2) +
  theme(legend.position = "top",
        axis.title = element_blank(),
        axis.line = element_line(colour = "grey50", size = 1),
        panel.background = element_blank(),
        legend.background = element_rect(fill = "lemonchiffon", 
                                         colour = "grey50", 
                                         size = 1),
        legend.title = element_blank())

## Women vs men

ddd_women_vs_men_by_period <- bind_rows(men_by_period %>% mutate(group="men"),
                women_by_period %>% mutate(group="women")) %>% 
  dplyr::select(time, group, ATT, CI.lower, CI.upper, se=S.E.) %>% 
  pivot_wider(., names_from=group, 
              values_from=c(ATT, CI.lower, CI.upper, se)) %>% 
  summarise(mean=round((ATT_women-ATT_men), 2),
            se=sqrt(se_men^2 + se_women^2),
            lower=round((mean-qnorm(p=1-(0.05)/2, mean=0, sd=1)*se), 2),
            upper=round((mean+qnorm(p=1-(0.05)/2, mean=0, sd=1)*se), 2),
            pvalue=round(pnorm(0,mean=abs(mean),sd=se) + 
              (1 - pnorm(0,mean=-abs(mean),sd=se)),2),
            time=time) %>% 
  relocate(time)

s2 <- ddd_women_vs_men_by_period %>% 
  ggplot(aes(x = time, y = mean, color = '#4daf4a')) +
  scale_color_manual(values=c('#4daf4a'), labels=c('Women vs Men'))+
  geom_point(position = position_dodge(width = .5)) + 
  geom_errorbar( aes( ymax = upper, ymin = lower) , width = .2, position = position_dodge(width = .5), size = 0.7 )+
  geom_line(position = position_dodge(width = .5), size = 0.7) +
  ylim(-38, 38)+
  geom_hline(yintercept = 0, lty=2) +
  scale_x_continuous(breaks=seq(0,6,1))+
  theme(legend.position = "top",
        axis.line = element_line(colour = "grey50", size = 1),
        panel.background = element_blank(),
        legend.background = element_rect(fill = "lemonchiffon", 
                                         colour = "grey50", 
                                         size = 1),
        legend.title = element_blank(),
        axis.title = element_blank())


triple_ddd <- ggarrange(s1,s2, ncol = 2, nrow=1) %>% 
  annotate_figure(left = textGrob("Change in CVD mortality, per 100,000 population", 
                                  rot = 90, vjust = 1, gp = gpar(cex = 1.3)), #cex = 0.8
                  bottom = textGrob("Number of years since Medicaid Expansion", gp = gpar(cex = 1.3)))

triple_ddd